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ABSTRACT 


-/ 

A  phase  convolution  filter  has  been  developed  for  the 
University  of  Washington's  ocean  bottom  seismometers  <OBS> 
to  correct  severe  waveform  distortion  on  playback.  Signals 
are  recorded  directly  on  magnetic  tape  after  compression  by 
taking  the  bipolar  square-root  of  the  signal.  Prior  to 
digitization!  the  analog  playback  head  shifts  the  phases  of 
the  «low  frequency  harmonics  and  thus  distorts  the  recorded 
waveform.  If  the  distorted  signal  is  squared  without 
correction*  the  original  waveforms  and  bandwidths  are  not 
recovered*  leading  to  errors  in  the  selection  of  arrival 
times  and  an  inability  to  analyse  the  seismic  coda  for  late 
arrivals.  The  phase-correcting  filter  has  been  constructed 
by  inverting  a  smoothed  version  of  the  spike  spectrum  that 
results  when  a  square  wave  is  recorded  on  the  tape. 
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Ocean  bottom  seismometers  have  been  used  in  marine 
geophysical  investigations  to  eliminate  the  shortcomings  of 
free  floating  surface  receivers.  The  University  of 
Washington's  OBS  as  described  by  Lister  and  Lewis*  1976.  was 
designed  with  simplicity  as  the  primary  goal  and  has  had  a 
very  high  recovery  rate  both  for  the  instruments  and  for 
usable  data  (Johnson  et  al.  1977).  We  will  describe  here 
only  those  features  of  the  OBS  pertinent  to  the  solution  of 
the  waveform  fidelity  problem  that  arose  when  the 
direct-recorded  data  tapes  were  played  back  into  a 
digitizer. 


Three  orthogonal  components  of  motion  are  measured 
using  L  IS  B  Mark  Products  4. 5  hertz  geophones.  They  are 
leveled  in  a  boat  floating  on  silicone  oil  of  a  thickness 
chosen  to  place  the  high  pass  filter  corner  frequency  well 
below  0.5  hertz.  The  signal  from  the  sensors  is  amplified 
and  passed  through  two  three-pole  Bessel  filters  to  define  a 
2  to  100  hertz  passband  with  a  minimum  of  time-envelope 
distortion.  The  filtering  limits  the  signal  to  within  the 
band  pass  capabilities  of  the  tape  recorder.  Dynamic-range 
compression  is  accomplished  by  taking  the  square  root  of  the 
signal  deviation.  negative  values  being  handled  by  a 
separate  circuit.  and  the  results  recombined  with  sign 
preservation.  Square-rooting  generates  low  harmonic  and 
cross-modulation  distortion.  consistent  with  substantial 
compression,  and  does  not  require  a  calibration  level,  since 
the  transition  to  linear  response  can  be  made  low  enough  to 
vanish  in  the  background  noise  level.  The  harmonic  content 
of  a  square-rooted  sine  wave  is  about  26  percent,  and  is 
independent  of  amplitude.  The  signal  is  recorded  by  a  slow 
speed  direct  method  on  tape  moving  at  1mm  per  second, 
together  with  a  clock  signal  encoded  with  the  day.  hour, 
minute,  and  second.  On  playback,  a  mini-computer  is  used  to 
digitize  all  four  channels  and  decode  the  clock  in  software. 
The  analog  tapes  are  played  back  on  a  slightly  modified 
Hewlett-Packard  3960  series  tape  recorder.  at  a  speed 
approximately  25  times  the  recording  speed.  A  four-pole  20 
hertz  notch  filter  is  used  to  eliminate  the  noise  induced  by 
vibration  of  the  OBS  tape  drive  motor.  The  signal  is 
amplified  by  a  tracking  three-channel  buffer  amplifer  until 
the  largest  signal  peaks  approach  the  limits  of  the 
digitizer's  range.  The  sample  spacing  is  0.01  seconds  in 
OBS  time,  or  2500  samples  per  second  at  playback  speed.  so 
that  a  1250  hertz  low-pass  Butterworth  filter  could  be 
inserted  in  the  analog  signal  processing  chain  to  reduce 
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aliasing. 


In  figure  1,  tha  measured  frequency  response  of  the 
electronic  package  (filled  circlas)  is  compared  to  the 
calculated  response  (open  circles)  of  the  electronics  plus 
geophone  sensor  (Lee  Bond/  unpublished  technical  report).  A 
2  volt  signal  was  inserted  at  the  calibration  input/ 
attenuated  by  1000/  and  then  measured  at  the  record  head  to 
produce  the  solid  line  in  figure  1.  The  theoretical 
response  of  the  geophone  was  calculated  from  the 
manufacturer 's  specifications  and  combined  with  the 
calculated  response  of  the  electronics.  The  combined 
response  is  minus  3db  in  amplitude  at  4.5  hertz/  but  is 
essentially  flat  from  8  tc  50  hertz.  The  square  root 
circuitry  was  calibrated  in  a  cold  room  to  verify  that  it 
follows  the  proper  power  law.  This  data  (figure  2)  produced 
the  regression  relationship  VOLTSout=SQRT(K*VOLTSin >  for  an 
input  voltage  of  10  millivolts  to  10  volts.  The  constant 
factor/  Kz  is  4.  5. 


A  5  hertz  calibration  signal  was  recorded  on  the  OBS/ 
digitized/  and  squared  to  show  the  type  of  waveform 
distortion  generated  by  the  analog  playback  head  (figure  3). 
In  addition/  the  digitized  output  of  a  perfect  2  hertz 
square  wave  (figure  4)  demonstrates  that  the  lowest  harmonic 
has  been  phase  shifted  but  the  higher  harmonics  which 
combine  to  preserve  the  sharp  onset/  have  not. 


FILTER  DESIGN 


A  filter  can  be  realized  by  comparison  of  a  distorted 
output  with  the  original  undistorted  signal/  in  frequency 
space/  simply  by  the  division  of  the  latter  by  the  former. 
The  frequency  domain  filter  can  then  be  inverted  to  a  time 
domain  filter/  and  applied  to  the  signal  by  convolution  to 
avoid  the  windowing  problems  of  finite  discrete  transforms. 
When  constructing  this  type  of  filter/  it  is  helpful  to  have 
the  discrete  frequencies  of  the  Fourier  transform  matched  to 
the  harmonics  present  in  the  test  signal.  This  criterion 
avoids  the  need  for  vectorial  addition  of  the  Fourier 
frequencies  near  the  test  wave  frequencies/  and  is  satisfied 
by  having  a  whole  number  of  test  wave  cycles  in  the 
transform  time  window. 
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FIGURE  1.  Measured  frequency  response  of  the  DBS.  The 
•  solid  line  is  the  response  of  the  electronics.  The  dashed 

line  is  the  calculated  response  of  the  electronics  including 
the  geophone  sensor.  Redrawn  from  L.  Bond.  unpublished 
technical  report  1979. 
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FIGURE  2.  Measured  transfer  function  of  the  OBS  electronics 
including  square-root  circuitry.  The  data  was  taken  at  5 
herti.  Redrawn  from  L.  Bond#  unpublished  technical  report 
1979 
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PROPERTIES  OF  FOURIER  TRANSFORM  PAIRS 


The  Fourier  pair*  f(t>  amd  F(w)  have  properties  that 
simplify  computation  in  tha  frequency  domain.  A  physically 
realizable  tima  function,  i. a.  any  gaophysical  record.  is  a 
pura  raal  function,  and  the  corresponding  F(w)  is  complex 
and  is  defined  for  both  positive  and  negative  w  up  to  the 
Nyquist  frequency.  Thus  it  contains  redundant  information, 
described  by 

This  implies  that  the  real  part  of  F(w)  is  even  (the  same 
for  positive  and  negative  w>  and  the  imaginary  part  is  odd. 
A  function  fulfilling  these  conditions  is  called  Hermitian. 
All  geophysical  spectra  are  Hermitian.  Furthermore,  since 

|FCufc\e^Cu* 

then  the  modulus  F(w)  is  even  and  the  phase  ^(w>  is  odd. 

The  time-shifting  theorem  is  also  a  useful  property. 
Given  that  f(t)  transforms  to  F(w>.  then  for  a  real  constant 
. 

If  a  time  function  is  shifted  by  a  constant  amount,  then  its 
amplitude  spectrum  remains  unchanged  while  its  phase 
spectrum  changes  by  a  linear  amount  determined  by  the  time 
shift.  This  property  is  useful  to  control  the  position  of 
the  final  time  domain  filter  in  the  transform  window. 
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EXPERIMENTAL  PROCEDURE 


An  experiment  was  designed  that  would  test  the  phase 
response  of  the  playback  electronics  within  the  seismic 
bandwidth  2  to  25  hertz.  The  test  signal  was  the  clock's 
internal  2  kilohertz  oscillator  that  was  divided  down  by 
1024  to  produce  an  accurate!  stable  square  wave  of  frequency 
1.9531250  hertz.  The  number  of  sample  points  used  in  the 

This  assured 


Fast  Fourier  Transform  (FFT)  routine  was  1024. 
a  whole  number  of  cycles>  20,  in  the  sample 
frequencies  present  in  a  square  wave  are  the 
plus  an  exponentially  decreasing  series 
harmonics.  Thus,  the  frequencies  sampled  using 
wave  are  approximately  2, 6, 10. .. hertz.  The  test  signal  was 
inserted  at  the  record  head  amplifier  of  all  three  channels 
but  by-passed  the  other  active  electronics  of  the  OBS.  The 
clock  signal  was  recorded  to  synchronize  the  digitization  on 
p laybac  k . 


window.  The 
fundamental 
of  the  odd 
this  square 


COMPUTATION  OF  THE  FILTER 


The  filter  is  constructed  in  the  frequency  domain 
utilizing  only  the  fundamental  and  odd  harmonic  frequencies. 
The  algebra  of  the  filter  is  as  follows.  One  can  define  the 
principal  time  functions  as; 

s(t)=input  square  wave 
d ( t )=d istorted  output 
f < t )«invtrs»  filter 
g(t)“forward  filter, 

and  in  frequency  space  each  time  function  Fourier  transforms 
to, 

s(t)  becomes  S(w) 
d(t)  becomes  D(w) 
f<t)  becomes  F(w). 


Simple  filter  theory  states  that  an  input  is  convolved  (< 
with  a  filter  to  produce  an  output  (equation  1). 


s ( t )*g  <  t )=d ( t  > 


(1) 
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An  invtrsi  filter,  f(t)»  can  ba  found  such  that 

g(t)*f (t)*delta(t>. 

the  unit  spike  function.  Convolving  f(t)  with  equation  1 
produces 

s( t)*g(t)*f (t)»d(t)#f (t). 

Since  a  function  convolved  with  a  delta  function  is  the 
original  function<  then 


s ( t  >*d ( t )*f ( t ) . 


( 2  > 


If  a  suitable  filter  can  be  constructed!  then  the  original 
signal  is  recovered  bg  convolving  f(t)  with  the  distorted 
output. 
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Convolution  in  one  domain  is  equal  to  multiplication  in 
the  transform  domain.  Thus  equation  2  can  be  directly 
written  in  Fourier  space  as< 

S(w)=D(w).  F(w) 
or  (3) 

F( w)*S(w) /D(w) . 

The  functions  S(w) .  0(w)< and  F(w>  are  complex  (equation  4)  so 
that  the  filter  F(w)  can  be  expressed  as  a  complex  division 
(equation  5). 

-SCujVXs  +  'v'/s  . 

QlwVX.+W/o  '  ' 

= r  + 1 

L  J  L  xS  +  yb'  J 

The  frequencies  present  in  the  test  square  wave  are  the 
odd  harmonics  of  1.933251  hertz.  These  frequencies  are 
identical  to  the  discrete  Fourier  frequencies  calculated  for 
1024  points  bg  the  FFT  program.  Thus  only  the  complex 
values  at  the  harmonics  are  used  in  equation  5  to  calculate 
the  complex  values  of  the  filter. 
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The  filter  is  also  Hermitian  and  must  satisfy  the 
oddness/eveness  properties  in  a  smooth  fashion.  Figure  5 
displays  two  partial  plots  of  the  imaginary  values  of  the 
filter/  near  the  edges  of  the  frequency  window  where  the 
amplitude  of  the  imaginary  part  must  go  through  zero.  The 
circles  are  the  discrete  frequency  values  and  the  connecting 
line  was  calculated  using  a  cubic  spline  interpolation 
routine.  However#  two  alterations  were  made  in  generating 
the  continuous  version.  First/  a  Gaussian  curve  was  applied 
from  0  to  6  hertz  with  the  1.95  hertz  value  eliminated. 
This  produced  a  smooth  transition  to  zero  hertz  and 
prevented  the  filter  from  trying  to  compensate  for  the 
natural  low  frequency  rolloff  of  the  playback  head.  At  the 
folding  point  or  Nyquist  frequency#  the  oddness  requirement 
was  satisfied  by  eliminating  the  48.8  hertz  value  and  by  a 
similar  Gaussian  smoothing  from  45  to  50  hertz.  The  real 
part  of  the  filter  was  cubic  spline  interpolated  between 
harmonic  values  with  a  smooth  continuation  to  the  folding 
point  (figure  6). 


The  functions  S(w)  and  D(w)  can  be  alternatively 
expressed  as  Fourier  series  (equation  6). 

SW  ~  2.  ^ e-  e*') 

OM  -  Dw=  4  0*^  1 

The  expression  for  the  filter  in  equation  3  becomes# 

~F<o-  2L  (Ift el . 

The  response  of  the  filter  is  the  ratio  of  the  amplitude 
responses#  Sn  and  Dn>  and  the  difference  in  the  phase 

responses.  The  amplitudes  at  the  harmonics  are  practically 

equal  (figure  7).  and  the  ratio  (Sn/Dn)  is  approx imately 

unity.  In  other  words#  the  filter  does  not  alter  the 

Fourier  amplitudes  in  the  signal  significantly. 


The  phase  response  of  a  perfect  square  wave  is  a 
constant.  However#  an  initial  time  shift  will  result  in  a 
linear  slope  in  the  phases.  This  slope  is  seen  in  the  phase 
responses  of  Sw  and  Dw  (figure  8).  The  difference  in  these 
phases  (lower  plot)  is  the  response  of  the  filter  but  with  a 
different  slope.  The  deviation  from  this  slope  is  the  phase 
angle  correction  that  is  applied  to  the  signal.  The  phase 
distortion  occurs  in  the  seismic  bandwidth  2  to  15  hertz 
(figure  8)  and  is  most  severe  below  10  hertz. 
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FIGURE  9.  Plots  of  the  imaginary  part  of  tho  complex 
invert*  filter.  Ordinate  units  are  arbitrary.  See  text  for 
d iscutsion. 
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FIGURE  7.  A  comparison  of  the  amplitude  spectra  of  a 
perfect  square  wave  (filled  circles)  with  the  spectra  of  the 
distorted  output  (crosses).  The  amplitude  scale  is  linear 
and  the  plots  have  been  normalised  to  1. 
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TIME  DOMAIN  REPRESENTATION  OF  THE  FILTER 


The  phase-correcting  filter  was  constructed  and 
modified  in  the  frequency  domain.  It  is  inverted  to  a  time 
domain  function  by  an  inverse  Fourier  transform  and  consists 
of  1024  paints,  the  same  number  as  in  the  frequency  domain. 
Figure  9  is  a  plot  of  the  middle  2S0  points  of  the  filter. 
The  energy  is  limited  to  the  center  points  with  the  tails 
being  near  zero.  The  filter's  amplitude  spectrum  can  be 
regarded  as  flat,  even  on  this  exaggerated  linear  scale 
(figure  9).  The  center  SO  points  (0.5  seconds)  of  the 
filter  retains  greater  than  99  percent  of  the  total  power  of 
the  original  filter.  This  length  filter  can  be  effectively 
convolved  with  the  signal  to  produce  the  desired  results 
In  addition.  the  resultant  SO  point  operator  is 
computational ly  efficient  for  moderate  length  data  sets. 


The  final  truncated  filter  is  smoothed  by  a  running 
3-point  binomial  filter.  This  has  the  desired  effect  of 
low-pass  filtering  the  inverse  filter  to  remove  any  residual 
high  frequency  oscillations.  However,  the  filter  remains 
flat  over  the  seismic  bandwidth  5  to  20  hertz.  and  the 
response  is  minus  3db  in  amplitude  at  3  and  22  hertz  with  a 
slow  roll-off  to  the  Nyquist  frequency.  50  hertz. 


RESULTS 


A  test  of  the  filter  was  made  by  convolving  the 
truncated  version  with  the  original  distorted  square  wave. 
Figure  10  is  a  sketch  of  several  cycles  of  the  filtered 
square  wave.  The  phase  correction  has  been  accomplished 
effectively.  A  very  slight  rounding  and  overshoot  is 
present  on  the  square  wave  edges  due  to  the  steps  taken  to 
smooth  the  filter,  but  the  anomalies  are  smaller  than  the 
random  variations  due  to  tape  and  head-ampl if ier  noise. 


The  final  result  obtained  by  convolving  the  phase 
correction  filter  with  a  seismic  trace  is  displayed  in 
figure  11  (trace  PHA8E).  The  initial  trace,  labeled  RAW.  is 
the  digitized.  square-rooted  signal.  The  SOR-trace  is  the 
result  of  bi-polar  squaring  the  RAW-signal.  The  waveform  is 
overall  spiky  in  nature  and  the  first-energy  arrival  time  is 
a  half-cycle  late,  about  0.05  seconds.  Additional  traces 
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FIGURE  9  A  linear  plot  of  the  amplitude  spectra  of  the 
inverse  phase  filter 
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FIGURE  10.  A  sketch  of  the  waveform  that  results  from 
convolving  the  distorted  square  wave  signal  with  the  inverse 
phase  filter. 
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FIGURE  11.  A  comparison  between  the  RAW  seismic  signal  and 
tha  phase  corractad  signal)  trace  PHASE.  The  SQR-trace  is 
the  raw  signal  squared,  and  tha  LP-tracas  are  various  low- 
pass  filtered  signals.  Refer  to  text  for  discussion. 
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were  generated  front  the  raw  squared  data  by  progressively 
filtering  out  more  and  more  of  the  high  frequencies/ 
consisting  mostly<  but  not  exclusively/  of  spurious 
harmonics  of  the  basic  seismic  signal.  They  are  designated 
LP-  and  the  -3db  cutoff  point  in  hertz  of  a  zero-phase  shift 
low-pass  convolution  filter  less  than  0.5  seconds  in  length. 
The  difference  between  these  LP~signals  and  trace-PHASE  are 
subtle  but  significant.  The  general  waveform  progresses 
from  being  very  spiky  (LP-25)  to  being  overly  'filtered' 
looking  (LP-10).  The  first  motion  of  the  PHASE-trace  is  a 
distinct  positive  (upward)  pulse.  Only  trace  LP-10  begins 
to  show  the  first  arrival  properly/  but  it  is  still  several 
hundreths  of  a  second  late.  It  should  be  noted  that  the 
amplitudes  for  each  trace  are  relative  to  a  maximum 
excursion  of  one  inch/  and  therefore  can  not  be  used  for 
direct  comparison.  However/  the  amplitude  of  the  first 
arrival  wave  packet  of  the  LP-traces  were  usually  no  more 
than  80  percent  of  that  of  the  PHASE-trace. 


CONCLUSIONS 


Until  recently  nearly  all  the  results  obtained  from 
marine  seismic  refraction  have  been  based  on  the  analysis  of 
the  travel-time  distance  relations  for  the  most  prominent 
arrivals  on  the  records.  In  this  report/  and  its  companion 
(  Determination  of  the  orientation  and  gain  of  horizontal 
geophones  used  in  ocean  bottom  seismometers/  Technical 
Report  Number  386  )  new  methods  have  been  described  which 
attempt  to  obtain  a  more  detailed  model  of  the  velocity 
structure  by  making  use  of  more  of  the  information  on  the 
seismograms.  The  combinations  of  travel-time  analysis  with 
studies  of  the  true  amplitudes/  particle  motions  and 
waveforms  of  the  arrivals/  places  much  greater  constraints 
on  the  velocity  structure  and  increases  the  potential 
resolution  of  the  refraction  technique.  In  the  future/  the 
spectral  content  of  an  arrival  may  be  as  valuable  in 
deducing  detailed  structure  as  the  timing  of  that  arrival/ 
or  its  raw  amplitude.  To  obtain  this  information/  it  is 
essential  that  the  characteristics  of  the  receiver  and 
playback  electronics  be  well  known  and  corrected  for  in 
processing  as  well  as  practicable. 
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The  unprocessed  signal  from  the  University  of 
Washington  OBS  showed  that  some  form  of  distortion  was 
occuring  in  the  seismic  bandwidth.  By  recording  a 
synchronised  square  wave,  and  comparing  the  output  with  the 
\nput  in  the  frequency  domain,  it  became  obvious  that  the 
lower  harmonics  of  the  square  wave  were  being  phase  shifted. 
A  phase  correcting  filter  was  constructed  in  the  frequency 
domain  and  inverted  to  a  time  domain  convolution  operator  of 
50  points  <0.5  second)  length. 
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